Langerhans

Data from isolated Langerhans islets with Chromium Next GEM Single Cell 5p RNA library v2. Aligned with STARsolo.

Load files

library(Seurat)
library(Matrix)

fpath <- '/home/vqf/proyectos/TFGs/Ruben/NoteBooks/mat/alignments/'
fsuffix <- 'Solo.out/Gene/filtered'
smpls <- c('158WT', '238WT', '244WT', '161KO', '235KO', '243KO')
types <- c('WT', 'WT', 'WT', 'KO', 'KO', 'KO')
tmpobjs <- c()
for (i in (1:length(smpls))){
  smpl <- smpls[i]
  tpe <- types[i]
  mtx <- paste(fpath, smpl, '/', smpl, fsuffix, '/matrix.mtx', sep = '')
  cll <- paste(fpath, smpl, '/', smpl, fsuffix, '/barcodes.tsv', sep = '')
  ftr <- paste(fpath, smpl, '/', smpl, fsuffix, '/features.translated.tsv', sep = '')
  interm <- ReadMtx(mtx = mtx, cells = cll, features = ftr)
  langer <- CreateSeuratObject(interm, project = 'Langerhans')
  langer$orig.ident <- smpl
  langer@meta.data$genotype <- tpe
  tmpobjs <- c(tmpobjs, langer)
  rm(interm)
}

langer <- merge(x = tmpobjs[[1]], y = tmpobjs[2:length(tmpobjs)], add.cell.ids=smpls)
str(langer)

Formal class 'Seurat' [package "SeuratObject"] with 13 slots
  ..@ assays      :List of 1
  .. ..$ RNA:Formal class 'Assay5' [package "SeuratObject"] with 8 slots
  .. .. .. ..@ layers    :List of 6
  .. .. .. .. ..$ counts.1:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. .. ..@ i       : int [1:12625864] 132 133 195 280 396 407 542 549 625 648 ...
  .. .. .. .. .. .. ..@ p       : int [1:9136] 0 1306 5962 7093 9605 10537 11822 13335 14320 16832 ...
  .. .. .. .. .. .. ..@ Dim     : int [1:2] 56980 9135
  .. .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. ..@ x       : num [1:12625864] 1 1 1 1 1 1 1 1 2 3 ...
  .. .. .. .. .. .. ..@ factors : list()
  .. .. .. .. ..$ counts.2:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. .. ..@ i       : int [1:12764382] 450 620 1088 1270 1533 1629 1656 2337 2492 2730 ...
  .. .. .. .. .. .. ..@ p       : int [1:9347] 0 209 970 1626 2549 3625 5719 6952 8045 8578 ...
  .. .. .. .. .. .. ..@ Dim     : int [1:2] 56980 9346
  .. .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. ..@ x       : num [1:12764382] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. .. .. .. .. ..@ factors : list()
  .. .. .. .. ..$ counts.3:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. .. ..@ i       : int [1:8332109] 149 239 251 455 495 556 558 625 628 712 ...
  .. .. .. .. .. .. ..@ p       : int [1:9225] 0 1187 2246 3194 3866 4500 5387 5856 6938 7669 ...
  .. .. .. .. .. .. ..@ Dim     : int [1:2] 56980 9224
  .. .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. ..@ x       : num [1:8332109] 1 1 1 2 1 1 1 1 1 5 ...
  .. .. .. .. .. .. ..@ factors : list()
  .. .. .. .. ..$ counts.4:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. .. ..@ i       : int [1:21186415] 125 185 239 251 271 319 481 490 503 556 ...
  .. .. .. .. .. .. ..@ p       : int [1:12284] 0 1235 2980 5220 6477 8583 9381 11054 12627 12875 ...
  .. .. .. .. .. .. ..@ Dim     : int [1:2] 56980 12283
  .. .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. ..@ x       : num [1:21186415] 1 1 1 4 1 1 2 2 1 3 ...
  .. .. .. .. .. .. ..@ factors : list()
  .. .. .. .. ..$ counts.5:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. .. ..@ i       : int [1:22457424] 87 149 150 172 233 239 251 271 274 278 ...
  .. .. .. .. .. .. ..@ p       : int [1:11602] 0 1566 2561 4913 5904 8185 9762 11143 12669 19210 ...
  .. .. .. .. .. .. ..@ Dim     : int [1:2] 56980 11601
  .. .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. ..@ x       : num [1:22457424] 1 1 1 1 2 1 3 7 1 1 ...
  .. .. .. .. .. .. ..@ factors : list()
  .. .. .. .. ..$ counts.6:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. .. ..@ i       : int [1:22697606] 22 149 206 228 229 251 301 322 380 427 ...
  .. .. .. .. .. .. ..@ p       : int [1:10124] 0 3062 5843 7576 9831 15634 21951 24022 27154 29339 ...
  .. .. .. .. .. .. ..@ Dim     : int [1:2] 56980 10123
  .. .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. .. ..$ : NULL
  .. .. .. .. .. .. ..@ x       : num [1:22697606] 1 1 1 1 1 2 1 3 2 2 ...
  .. .. .. .. .. .. ..@ factors : list()
  .. .. .. ..@ cells     :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
  .. .. .. .. .. ..@ .Data: logi [1:61712, 1:6] TRUE TRUE TRUE TRUE TRUE TRUE ...
  .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. ..$ : chr [1:61712] "158WT_AAACCTGAGTGGTCCC" "158WT_AAACCTGCAATGCCAT" "158WT_AAACCTGCACCCATGG" "158WT_AAACCTGCAGATCTGT" ...
  .. .. .. .. .. .. .. ..$ : chr [1:6] "counts.1" "counts.2" "counts.3" "counts.4" ...
  .. .. .. .. .. ..$ dim     : int [1:2] 61712 6
  .. .. .. .. .. ..$ dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:61712] "158WT_AAACCTGAGTGGTCCC" "158WT_AAACCTGCAATGCCAT" "158WT_AAACCTGCACCCATGG" "158WT_AAACCTGCAGATCTGT" ...
  .. .. .. .. .. .. ..$ : chr [1:6] "counts.1" "counts.2" "counts.3" "counts.4" ...
  .. .. .. ..@ features  :Formal class 'LogMap' [package "SeuratObject"] with 1 slot
  .. .. .. .. .. ..@ .Data: logi [1:56980, 1:6] TRUE TRUE TRUE TRUE TRUE TRUE ...
  .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. ..$ : chr [1:56980] "Gm37671" "Gm19087" "Gm8941" "Gm38212" ...
  .. .. .. .. .. .. .. ..$ : chr [1:6] "counts.1" "counts.2" "counts.3" "counts.4" ...
  .. .. .. .. .. ..$ dim     : int [1:2] 56980 6
  .. .. .. .. .. ..$ dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:56980] "Gm37671" "Gm19087" "Gm8941" "Gm38212" ...
  .. .. .. .. .. .. ..$ : chr [1:6] "counts.1" "counts.2" "counts.3" "counts.4" ...
  .. .. .. ..@ default   : int 6
  .. .. .. ..@ assay.orig: chr(0) 
  .. .. .. ..@ meta.data :'data.frame': 56980 obs. of  0 variables
  .. .. .. ..@ misc      : list()
  .. .. .. ..@ key       : chr "rna_"
  ..@ meta.data   :'data.frame':    61712 obs. of  4 variables:
  .. ..$ orig.ident  : chr [1:61712] "158WT" "158WT" "158WT" "158WT" ...
  .. ..$ nCount_RNA  : num [1:61712] 3918 24759 2859 7331 1551 ...
  .. ..$ nFeature_RNA: int [1:61712] 1306 4656 1131 2512 932 1285 1513 985 2512 912 ...
  .. ..$ genotype    : chr [1:61712] "WT" "WT" "WT" "WT" ...
  ..@ active.assay: chr "RNA"
  ..@ active.ident: Factor w/ 1 level "Langerhans": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "names")= chr [1:61712] "158WT_AAACCTGAGTGGTCCC" "158WT_AAACCTGCAATGCCAT" "158WT_AAACCTGCACCCATGG" "158WT_AAACCTGCAGATCTGT" ...
  ..@ graphs      : list()
  ..@ neighbors   : list()
  ..@ reductions  : list()
  ..@ images      : list()
  ..@ project.name: chr "SeuratProject"
  ..@ misc        : list()
  ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
  .. ..$ : int [1:3] 5 1 0
  ..@ commands    : list()
  ..@ tools       : list()

rm(tmpobjs)

counts1 <- GetAssayData(object = langer, layer = "counts.1")
counts1[200:210, 1:3]

11 x 3 sparse Matrix of class "dgCMatrix"
          158WT_AAACCTGAGTGGTCCC 158WT_AAACCTGCAATGCCAT 158WT_AAACCTGCACCCATGG
Epha4                          .                      .                      .
Mir6352                        .                      .                      .
Colgalt2                       .                      .                      .
Gm28386                        .                      .                      .
Gm28385                        .                      .                      .
Gm26263                        .                      .                      .
Gm28387                        .                      .                      .
Poglut2                        .                      .                      .
Btf3-ps14                      .                      .                      .
Catspere2                      .                      .                      .
Gm37313                        .                      .                      .

Now all objects are in langer. Starting sample is in langer@meta.data$orig.ident and genotype is in langer@meta.data$genotype.

Pre-processing

table(langer@meta.data$genotype)

   KO    WT 
34007 27705 

Mitochondrial

langer <- PercentageFeatureSet(langer, pattern = "^mt-", col.name = "percent.mt")
table(langer@meta.data$genotype)

   KO    WT 
34007 27705 

Visualizations and filtering

Visualize channels for claim

VlnPlot(langer, features = c("nFeature_RNA"), group.by = "orig.ident", pt.size = 0, y.max = 7500)

VlnPlot(langer, features = c("nCount_RNA"), group.by = "orig.ident", pt.size = 0, y.max = 7500)

VlnPlot(langer, features = c("percent.mt"), group.by = "orig.ident", pt.size = 0)

VlnPlot(langer, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), group.by = "orig.ident", ncol = 3, pt.size = 0)

VlnPlot(langer, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), group.by = "genotype", ncol = 3, pt.size = 0)

Genes per cell

reads_per_cell <- langer[['nFeature_RNA']]$nFeature_RNA
MIN_READS_PER_CELL <- 300
MAX_READS_PER_CELL <- 6000
MAX_MT <- 5
plot(sort(reads_per_cell), xlab='cell', log='y', main='Reads per cell (ordered)')
abline(h=MIN_READS_PER_CELL, col='magenta')

abline(h=MAX_READS_PER_CELL, col='blue')

table(langer@meta.data$genotype)

   KO    WT 
34007 27705 

langer <- subset(langer, percent.mt<MAX_MT)
langer.highexpr <- subset(langer, nFeature_RNA>MAX_READS_PER_CELL)
langer <- subset(langer, nFeature_RNA>MIN_READS_PER_CELL)
langer <- subset(langer, nFeature_RNA<MAX_READS_PER_CELL)
VlnPlot(langer, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), group.by = "orig.ident", ncol = 3, pt.size = 0)

VlnPlot(langer, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), group.by = "genotype", ncol = 3, pt.size = 0)

Normalize and scale

langer <- SCTransform(langer, vars.to.regress = "percent.mt", verbose = FALSE)
table(langer@meta.data$genotype)

   KO    WT 
31588 26893 

langer <- RunPCA(langer) 
DimHeatmap(langer, dims = 1:10, cells = 50, balanced = TRUE)

ElbowPlot(langer)

langer <- RunTSNE(langer, dims = 1:20)
langer <- RunUMAP(langer, dims = 1:20)

Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

TSNEPlot(langer, split.by="genotype")

FeaturePlot(langer, features = c("Gcg", "Ins1", "Ins2", "Jchain", "Sst", "Ppy", "Ghrl", "Cd68", "Cd4", "Pecam1", "Tagln", "Iapp", "Col3a1"), cols = c("grey", "red"), reduction = "tsne", combine = F, split.by = "orig.ident")

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

[[19]]

[[20]]

[[21]]

[[22]]

[[23]]

[[24]]

[[25]]

[[26]]

[[27]]

[[28]]

[[29]]

[[30]]

[[31]]

[[32]]

[[33]]

[[34]]

[[35]]

[[36]]

[[37]]

[[38]]

[[39]]

[[40]]

[[41]]

[[42]]

[[43]]

[[44]]

[[45]]

[[46]]

[[47]]

[[48]]

[[49]]

[[50]]

[[51]]

[[52]]

[[53]]

[[54]]

[[55]]

[[56]]

[[57]]

[[58]]

[[59]]

[[60]]

[[61]]

[[62]]

[[63]]

[[64]]

[[65]]

[[66]]

[[67]]

[[68]]

[[69]]

[[70]]

[[71]]

[[72]]

[[73]]

[[74]]

[[75]]

[[76]]

[[77]]

[[78]]

VlnPlot(langer, features = c("Gcg", "Ins1", "Ins2", "Gapdh", "Ghrl"), group.by = "genotype", pt.size = 0) 

Get clusters and save object

langer <- FindNeighbors(langer, reduction = 'pca', dims = 1:20)
langer <- FindClusters(langer, resolution = 0.5)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 58481
Number of edges: 2037766

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9553
Number of communities: 28
Elapsed time: 12 seconds

DimPlot(langer, reduction = 'tsne', split.by =  'genotype', label = T, repel = T)

Cluster 17 probably contains beta and delta cells ###Get the genes used to classify the clusters

options(future.globals.maxSize= 8288490189)
langer <- PrepSCTFindMarkers(langer, assay = "SCT", verbose = FALSE)

all.markers <- FindAllMarkers(langer)
View(all.markers)
saveRDS(langer, "/home/vqf/proyectos/TFGs/Ruben/NoteBooks/langer.rds")
saveRDS(all.markers, "/home/vqf/proyectos/TFGs/Ruben/NoteBooks/allmarkers.rds")